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Abstract. We use a standard Monte-Carlo algorithm to study the slow 
dynamics of a binary Lennard-Jones glass-forming mixture at low temperature. 
We find that Monte-Carlo is by far the most efficient way to simulate a 
stochastic dynamics since relaxation is about 10 times faster than in Brownian 
Dynamics and about 30 times faster than in Stochastic Dynamics. Moreover, 
the average dynamical behaviour of the system is in quantitative agreement with 
the one obtained using Newtonian dynamics, apart at very short times where 
thermal vibrations are suppressed. We show, however, that dynamic fluctuations 
quantified by four-point dynamic susceptibilities do retain a dependence on the 
microscopic dynamics, as recently predicted theoretically. 



PACS numbers: 64.70.Pf, 05.20.Jj 

1. Introduction 

Numerical simulations play a major role among studies of the glass transition since, 
in contrast to experiments, the individual motion of a large number of particles can 
be followed at all times pQ. With present day computers, it is possible to follow 
the dynamics of a simple glass-forming liquid over more than 8 decades of time, 
and over a temperature window in which average relaxation timescales increase by 
more than 5 decades. However, at the lowest temperatures studied, relaxation is still 
orders of magnitude faster than in experiments performed close to the glass transition 
temperature. Nevertheless, it is now possible to numerically access temperatures 
which are low enough that many features associated to the glass transition physics 
can be observed: Strong decoupling phenomena [21 CUE], clear deviations from fits to 
the mode-coupling theory [5] (which are experimentally known to hold only at high 
temperatures), and crossovers towards activated dynamics |H1[7|. 

Computer simulations usually study Newtonian Dynamics (ND) by solving 
a discretized version of Newton's equations for a given pair interaction between 
particles [B]- Here, we study a glass- forming model in which a binary mixture of small 
and large particles interact via a Lennard-Jones pair potential, a model introduced 
by Kob and Andersen (KA) 0. It can be interesting to study also different types 
of microscopic dynamics for the same pair potential. If dynamics satisfies detailed 
balance with respect to the Boltzmann distribution, all structural quantities remain 
unchanged, although the dynamics might be very different. In colloidal glasses, for 
instance, the particles undergo Brownian motion arising from collisions with the 
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molecules of the solvent, and a stochastic dynamics is more appropriate. Theoretical 
considerations also suggest the study of different dynamics. Gleim et al. studied a 
Stochastic Dynamics (SD) to investigate whether the relaxation of the KA binary 
mixture depend on its microscopic dynamics, their answer being "no" UJ. In SD, a 
friction term and a random noise are added to Newton's equations, the amplitude of 
both terms being related by a fluctuation-dissipation theorem. Szamel and Flenner 
recently used Brownian Dynamics (BD) to study the same KA mixture ^H]- In this 
description there are no momenta, and positions evolve with a Langevin dynamics. 
They again find that relaxation using BD is very similar to the one resulting from 
ND. They emphasize that even the deviations from mode-coupling fitting are similar 
in BD and ND and conclude that momenta play no role in avoiding the mode- 
coupling singularity, contrary to previous claims but in agreement with more 
recent ones [T2"] . 

Recently, it was also discovered that dynamic heterogeneity, that is, spatio- 
temporal fluctuations around the average dynamical behaviour, sensitively depends 
upon the microscopic dynamics [HI 171 IT5]. In particular, a major role is played by 
conservation laws for energy and density. In the case of energy the mechanism can 
be physically understood as follows. For a rearrangement to take place in the liquid, 
the system has to locally cross an energy barrier. If the dynamics conserves the 
energy, particles involved in the rearrangment must borrow energy to the neighboring 
particles. This 'cooperativity' might be unnecessary if energy can be locally supplied 
to the particles by an external heat bath. Conservation laws, therefore, might 
introduce dynamic correlations between particles and dynamic fluctuations can be 
different when changing from Newtonian energy conserving dynamics to a stochastic 
thermostatted dynamics. This predicted influence of the microscopic dynamics on 
dynamic fluctuations [HJ d was in fact our principal motivation for the present study. 

In this article, we propose a third type of stochastic dynamics for the KA mixture 
and study in detail the dynamics of the system subjected to a standard Monte-Carlo 
(MC) dynamics. We find that MC is particularly efficient at relaxing the system since 
it is about 10 times faster than BD and 30 times faster than SD, while the average 
dynamics is still in quantitative agreement with ND. We are therefore in position to 
study both the very low temperature average dynamics of the model and its dynamic 
fluctuations in detail, shedding new light on both aspects. 

The paper is organized as follows. In Section|21we give details about the simulation 
technique and compare its efficiency to previously studied dynamics. In Section we 
present our numerical results. Section 0] concludes the paper. 

2. An efficient simulation technique 

We study a binary Lennard- Jones mixture made of Na = 800 and Nb — 200 particles 
of types A and B, respectively. Particles interact with the following Lennard- Jones 
pair potential 



where a, (3 6 [A, B] and r is the distance between the interacting pair of particles. 
Interaction parameters e a p and o~ a p are chosen to prevent crystallization and can be 
found in Ref. py. The length and energy are given in the standard Lennard- Jones 
units o~ aa (particle diameter), and caa (interaction energy), where the subscript A 




(1) 
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Figure 1. Self-intermediate scattering function, Eq. J2Ji at T = 0.5 and k = 7.21 
for various values of <5 max . Inset: The evolution of the relaxation time with <5 max 
unambiguously defines an optimal value <5 max ~ 0.15 for our simulations. 



refers to the majority species. The potential is truncated and shifted at a distance 
r = 2.5. Previous work ^ [S] has shown that the dynamics becomes slow below 
T « 1.0, while the fitted mode-coupling temperature for this system is T c « 0.435, 
although deviations from mode-coupling behaviour become noticable already below 
T « 0.47. 

We have implemented a standard Monte-Carlo dynamics |S] for the pair potential 
in Eq. (JTJ. An elementary move can be described as follows. A particle, i, located at 
the position is chosen at random. The energy cost, AEi, to move particle i from 
position Vi to a new position + Sr is evaluated, Sr being a random vector comprised 
in a cube of linear length S max centered around the origin. The Metropolis acceptance 
rate, p = min(l, e~~ @ *), where (3 = 1/T is the inverse temperature, is then used 
to decide whether the move is accepted. In the following, one Monte-Carlo timestep 
represents N — Na + Nb attempts to make such an elementary move, and timescales 
are reported in this unit. 

The one degree of freedom that remains to be fixed is <5 max which determines the 
average lengthscale of elementary moves. If chosen too small, energy costs are very 
small and most of the moves are accepted, but the dynamics is very slow because 
it takes a long time for particles to explore their cage. On the other hand too large 
displacements will on average be very costly in energy and acceptance rates can become 
prohibitively small. We seek a compromise between these two extremes by monitoring 
the dynamics at a moderately low temperature, T — 0.5, for several values of 5 max . As 
the most sensitive indicator of the relaxational behaviour we measure the contribution 
from the majority specie A to the self-intermediate scattering function, 

Fs (k,t) = ^p^w-^y 

We spherically average over wavectors of comparable magnitude, and present results 
for |k| = 7.21, which corresponds to the first diffraction peak in the static structure 
factor of the liquid. In Fig. U we present our results for £ max values between 0.05 and 
0.4. As expected we find that relaxation is slow both at small and large values of <5 max , 
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and most efficient for intermediate values. Interestingly we also note that the overall 
shape of the self-intermediate scattering function does not sensitively depend on <5 max . 

We define a typical relaxation time as F s (k,T a ) = e _1 and show its <5 max 
dependence in the inset of Fig. ^ A clear minimum is observed at the optimal value 
of <5max ~ 0.15. In the rest of the paper we only present data obtained for this value. 

As compared to previously studied dynamics, we find that, when expressed in 
numbers of integration timesteps, structural relaxation in Monte-Carlo simulations is 
marginally faster than in Newtonian dynamics, but 30 times faster than in Stochastic 
Dynamics 9 , and 10 times faster than in Brownian Dynamics ^01- We conclude 
therefore that MC is by far the most efficient way to perform stochastic molecular 
simulations of the present glass-forming material. 

The relative inefficiency of both BD and SD is due to the stochastic nature of their 
microscopic equations of motion. It is well-known that small integration timesteps are 
required for accurate integration of stochastic equations of motion, in particular to 
maintain the delicate balance between friction and noise required for the system to 
converge towards the correct equilibrium distribution |H] . No such constraint exists for 
MC dynamics, where elementary moves can be made arbitrarily large. Equilibrium 
only requires detailed balance to be fulfilled, and this is always the case with the 
Metropolis algorithm described above. With larger elementary moves, particles can 
efficiently explore their cage and relaxation is much faster. This physical interpretation 
is also supported by the optimal value <5 max = 0.15 that we report, which corresponds 
to a mean-square displacement of 0.225, very close to the plateau observed in the 
mean-square displacement shown in Fig. El (see below), which can be taken as a 
rough estimate of the cage size. Monte-Carlo simulations can of course be made 
even more efficient by implementing for instance swaps between particles, or using 
parallel tempering. The dynamical behaviour, however, is then strongly affected by 
such non-physical moves and only equilibrium thermodynamics can be studied. Since 
we want to conserve a physically realistic dynamics, we cannot use such improved 
schemes. 

We have performed simulations at temperatures between T = 2.0 and T = 0.43, 
the latter being smaller than the fitted mode-coupling temperature. For each 
temperature we have simulated 10 independent samples to improve the statistics. 
Initial configurations were taken as the final configurations obtained from previous 
work performed with ND [HE], so that production runs could be started immediately. 
For each sample, production runs lasted at least 15t q (at T — 0.43), much longer for 
higher temperatures. 

3. Results 

3.1. Average dynamics 

The self-intermediate scattering function, Eq. @, is shown in Fig. |2Jfor temperatures 
decreasing from T = 2.0 down to T = 0.43. These curves present well-known 
features. Dynamics at high temperature is fast and has an exponential nature. When 
temperature is decreased below T « 1.0, a two-step decay, the slower being strongly 
non-exponential, becomes apparent. Upon decreasing the temperature further, the 
slow process dramatically slows down by about 5 decades, while clearly conserving 
an almost temperature-independent non-exponential shape, as already reported for 
ND 0. 
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Figure 2. Left: Self-intermediate scattering function, Eq. J5Ji f° r = 7-21 and 
temperatures T = 2.0, 1.0, 0.75, 0.6, 0.5, 0.47, 0.45, 0.435, and 0.43 (from left to 
right). Right: Mean-squared displacement, Eq. J3J, for the same temperatures in 
the same order. 



Finally, as reported for SD 9 , we find that also the first process, the decay 
towards a plateau, slows down considerably when decreasing temperature. This 
process, called 'critical decay' in the language of mode-coupling theory is n °t 
observed when using ND, because it is obscured by the thermal vibrations occuring 
at high frequencies. Although the plateau seen in F s (k,t) is commonly interpreted 
as 'vibrations of a particle within a cage', the data in Fig. |3 discard this view. From 
direct visualisation of the particles' individual dynamics it is obvious that vibrations 
take place in just a few MC timesteps, while the decay towards the plateau can be as 
long as 10 4 time units at the lowest temperatures studied here. This decay is therefore 
necessarily more complex, most probably cooperative in nature. This interpretation 
is supported by recent theoretical studies where a plateau is observed in two-time 
correlators of lattice models where local vibrations are indeed completely absent ^5] . 
A detailed atomistic description of this process has not yet been reported, but would 
indeed be very interesting. 

Next, we study the mean-squared displacement for the majority specie. It is 
defined as 

1 n a 

i—l 

and we present its temperature evolution in Fig. [21 which mirrors the evolution of 
the self-intermediate scattering function in the same figure. Since we are studying 
a stochastic dynamics, displacements are diffusive at both short and long timescales. 
The plateau observed in F s (fc, t) now translates into a sub-diffusive regime in the mean- 
squared displacements separating the two diffusive regimes. At the lowest temperature 
studied, when t changes by three decades from 2 x 10 2 to 2 x 10 5 , the mean-squared 
displacement changes by a mere factor 2.2 from 0.02 to 0.044. Particles are therefore 
nearly arrested for several decades of times, before eventually entering the diffusing 
regime which allows for the relaxation of the structure of the liquid. 

3.2. Comparison to Newtonian and Stochastic Dynamics 

The previous subsection has shown that the Monte-Carlo dynamics of the KA mixture 
is qualitatively similar to the one reported for ND, apart at relatively short times where 
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Figure 3. Left: Comparison of the self-intermediate scattering function for 
k = 7.21, and T = 0.45, obtained in Monte-Carlo (MC) dynamics in this work, 
Newtonian Dynamics (ND) in Ref. and Stochastic Dynamics (SD) in Ref. 
Time is rescaled to obtain maximum overlap at large times. MC and SD agree 
over the complete time range (indeed the SD dotted line is barely visible below 
the full MC line), while MC and ND only agree when F B (k,t) is close to the 
plateau and below. The dip observed at short-time in ND is due to thermal 
vibrations suppressed in both SD and MC. Right: Temperature evolution of the 
alpha-relaxation time T a (T) and the inverse of the self-diffusion constant X/D(T) 
in an Arrhenius plot. Open symbols are for ND, closed symbols for MC (vertically 
shifted to obtain maximum overlap with ND data) the dashed lines are power law 
fits to a divergence at T c = 0.435, as originally reported in Ref. 



the effect of thermal vibrations is efficiently suppressed. We now compare our results 
more quantitatively with the dynamical behaviour observed using ND. 

In Fig. we compare the time dependence of the self-intermediate scattering 
function for three types of dynamics: the present Monte-Carlo data, the Newtonian 
Dynamics data taken from Ref. & , and the Stochastic Dynamics results from Ref. 9 , 
all obtained for the same parameters k = 7.21 and T = 0.45. We have rescaled the 
time to obtain maximum overlap in the long-time relaxation of the three curves. Quite 
strikingly, SD and MC data perfectly overlap over the complete time-range (8 decades 
of time) of the simulation. Indeed the SD dotted line is barely visible below the full 
line of the MC data in Fig. [3] This confirms our claim that MC defines a physically 
relevant microscopic dynamics, since it is completely equivalent to SD with the major 
advantage that it is 30 times faster, at least for the KA mixture. 

In Fig.0 we also confirm that the approach to the plateau is different in MC/SD 
and ND. In the latter, phonon-like vibrations affect the initial decay of F s (k,t). For 
instance, a shallow dip, generally attributed to the 'Boson peak', is observed at low 
temperature in ND, see the dashed line in Fig. [3] The long-time decay of the self- 
intermediate scattering function, however, is in full quantitative agreement for the 
three dynamics. This agreement was the main claim of Ref. [H], extended to BD in 
Ref. ^0] and for MC in the present work. 

Since all dynamics display similar long-time relaxation, it is sensible to also 
quantitatively compare the temperature evolution of the relaxation times, T a (T), 
already defined above. This is done in Fig. where we use a standard representation 
where an Arrhenius slowing down over a constant energy barrier, r Q ~ exp(i?/T), 
would appear as a straight line. The data clearly show some upwards bending in 
Fig. |3 which places the KA mixture in the family of fragile (though very weakly) glass- 
formers. We find that the temperature evolution of the alpha-relaxation time measured 
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in MC simulations is in complete quantitative agreement with the one obtained from 
ND, over the complete temperature range T — 2.0 — * 0.43. In particular the quality of 
a power-law fit of the slowing down, r a ~ (T — T c )~ 7 , as suggested by mode-coupling 
theory, is similar for both dynamics .5;, 9] . We have shown such a fit through our data, 
using the value T c = 0.435 determined in Ref. [5J- The fit describes the data over 
about 2.5 decades. Deviations from the mode-coupling fit appear below T w 0.47, and 
become obvious when T c is approached further. 

In Fig. |3| we also show the temperature evolution of the self-diffusion constant, 
defined from the long-time limit of the mean-square displacement as 

D=Um^. (4) 

The behaviour of the diffusion constant is qualitatively very close to the one of the 
alpha-relaxation time, and all the above remarks apply. The well-known difference 
between the two quantities is a slightly stronger temperature evolution of r a , implying 
a well-studied decoupling between translational diffusion and structural relaxation in 
this system 00], which is therefore very similar for different types of dynamics. 

Theoretically, an identical relaxation within MC/SD/BD/ND is an important 
prediction of mode-coupling theory ^3] because the theory uniquely predicts the 
dynamical behaviour from static density fluctuations. Gleim et al. argue that their 
finding of a quantitative agreement between SD and ND is a nice confirmation of this 
non-trivial mode-coupling prediction ;9< . Szamel and Flenner ^U] confirmed this claim 
using BD, and argued further that even deviations from mode-coupling predictions are 
identical. We confirm the validity of this statement even below T c , showing that the 
agreement between different dynamics, although indeed predicted by mode-coupling 
theory, is certainly valid at a much more general level. Similarly to Szamel and 
Flenner, we note that deviations from a power law divergence cannot be attributed 
to coupling to currents which are expressed in terms of particle velocities. In our MC 
simulations we have no velocities, so that avoiding the mode-coupling singularity is 
not due to the hydrodynamic effects pointed out in Ref. (see Ref. ^2] for more 
recent theoretical viewpoints). 



3.3. Multi-point susceptibility 

Having established the ability of MC simulations to efficiently reproduce the average 
slow dynamics obtained from ND simulations we now turn to the study of the dynamic 
fluctuations around the average dynamical behaviour, i.e. to dynamic heterogeneity. 

Dynamic fluctuations can be studied through the four-point susceptibility, x^if), 
which quantifies the strength of the spontaneous fluctuations around the average 
dynamics by their variance, 

xm = N A [(f*(Kt))-F?Qi,t)}, (5) 

where / s (k, t) = A^ 1 Y^j cos(k • [rj(t) — rj(0)]) represents the real part of the 
instantaneous value of the self- intermediate scattering function, so that F s (k,t) = 
(/ s (k, t)). As shown by Eq. (|5j), it is clear that Xi{t) wm be large if run-to-run 
fluctuations of the self-intermediate scattering functions are large. This is the case 
when the local dynamics becomes spatially correlated, as already discussed in several 

papers [HO El CHI US EHEH 

We show the time dependence of the dynamic susceptibility Xi {t) obtained from 
our MC simulations for various temperatures in Fig. 0] As predicted theoretically in 
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Figure 4. Four-point susceptibility, Eq. for the same temperatures as in 
Fig. |3 decreasing from left to right. We have highlighted with open circles the 
data at T = 0.43, which are fitted with two power laws shown as dashed lines 
with exponents 0.35 and 0.75 at short and large times, respectively. 

Ref. [23j we find that X4(£) presents a complex time evolution, closely related to the 
time evolution of the self-intermediate scattering function. Overall, Xi(t) ^ small at 
both small and large times when dynamic fluctuations are small. There is therefore 
a clear maximum observed for times comparable to r Q , where fluctuations are most 
prominent. The position of the maximum then shifts to larger times when temperature 
is decreased, tracking the alpha-relaxation. The most important physical information 
revealed by these curves is the fact that the amplitude of the peak grows when the 
temperature decreases. This is direct evidence that spatial correlations grow when the 
glass transition is approached. 

The two-step decay of the self-intermediate scattering function translates into a 
two-power law regime for Xi(t) approaching its maximum. We have fitted these power 
laws, Xi(t) ~ followed by Xi(f) ~ t h with the exponents a = 0.35 and b = 0.75 in 
Fig. 21 We have intentionally used the notation a and b for these exponents which are 
predicted, within mode-coupling theory, to be equal to the standard exponents also 
describing the time dependence of intermediate scattering functions Q3| . Our findings 
are in good agreement with previously reported values for a and b. See Refs. I2U] 
for a more extensive discussion and comparison to other theoretical predictions. 

We finally compare the dynamic susceptibility for various dynamics. In Fig. 
we present the time evolution of Xi(t) f° r a given temperature, T = 0.45 and four 
different dynamics: The present MC data, data from SD obtained in Ref. 6 , data 
for ND in the microcanonical (NVE) ensemble from Ref. jJJ, and data for ND in 
the canonical (NVT) ensemble from Ref. p]. To perform this comparison, we have 
again rescaled times to obtain the maximum overlap in the long-time region. In Fig. 03 
it is obvious that three curves are identical: ND-NVE, MC and SD data perfectly 
overlap near the maximum of Xi(t) an d have similar time dependences, apart at very 
short-times. On the other hand, ND-NVT data display a different time dependence 
and reveal considerably larger dynamic fluctuations in the long-time regime. 

We conclude therefore that, contrary to the average dynamics, the dynamic 
fluctuations quantified through the four-point susceptibility do retain a dependence 
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Figure 5. Four-point susceptibility for various dynamics and ensembles at 
T = 0.45. As in Fig. |3] times have been rescaled to obtain the maximum overlap 
in the long-time regime. The three overlapping thin lines represent data for ND- 
NVE dynamics, SD, and MC, while the thick line represents ND-ATVT data, 
for which dynamic fluctuations are clearly larger, as predicted theoretically and 
discussed in Ref. UJ. 



upon the microscopic dynamics since canonical estimates of Xi(f) are different for ND 
and for MC/SD/BD. Although perhaps counter-intuitive at first sight we find that 
dynamics with a stochastic heat-bath display dynamic fluctuations similar to the ones 
measured using microcanonical ND, while fluctuations are much larger in canonical ND 
simulations. As mentioned in the introduction, this confirms the idea that the energy 
conservation (implied by Newton's equations of motion) might lead to an amplification 
of dynamic fluctuations. With hindsight, this is not such a surprising result: The 
specific heat, after all, also behaves differently in different statistical ensembles. The 
ensemble dependence and dependence upon the microscopic dynamics are the main 
subjects of two recent papers jnj[7|. 

There is an experimentally relevant consequence of these findings. The difference 
between the microcanonical and canonical values of the dynamic fluctuations in ND 
can be shown to be equal to p|] 



where cy is the constant volume specific heat expressed in ks units. As shown in 
Fig. |S] the temperature derivative in Eq. JBJ represents in fact the major contribution 
to X4 VT ') meaning that the term Xa VE can be neglected in Eq. ©. Since the right 
hand side of © is more easily accessible in an experiment than Xi itself, Eq. © 
opens the possibility of an experimental estimate of the four-point susceptibility. This 
finding, and its experimental application to supercooled glycerol and hard sphere 
colloids, constitute the central results of Ref. |13j . 

4. Conclusion 




(6) 



We have implemented a standard Monte-Carlo dynamics on the well-known binary 
Lennard- Jones mixture introduced by KA. We have shown that the resulting average 
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dynamics is in full quantitative agreement with results from Newtonian dynamics, 
while being considerably faster than previously studied stochastic dynamics, namely 
Brownian and Stochastic dynamics. We have therefore at our disposal an efficient 
numerical technique to simulate the stochastic dynamics of the KA mixture at low 
temperature. This allowed us to show, in particular, that dynamic fluctuations retain a 
dependence upon the microscopic dynamics since four-point dynamical susceptibilities 
evaluated in the canonical ensemble for ND and MC quantitatively differ, because the 
energy conservation of Newton's equations amplify dynamic fluctuations. 
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